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Limiting eccentricity of sub-parsec massive black hole binaries 
surrounded by self-gravitating gas discs 
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ABSTRACT 

We study the dynamics of supermassive black hole binaries embedded in circumbinary 
gaseous discs, with the SPH code Gadget-2. The sub-parsec binary (of total mass M and 
mass ratio q = 1/3) has excavated a gap and transfers its angular momentum to the self- 
gravitating disc (Mdi S c = 0.2M). We explore the changes of the binary eccentricity e, by 
simulating a sequence of binary models that differ in the initial eccentricity eo, only. In ini- 
tially low-eccentric binaries, the eccentricity increases with time, while in high-eccentric bi- 
naries e declines, indicating the existence of a limiting eccentricity e cr i t that is found to fall 
in the interval [0.6, 0.8]. We also present an analytical interpretation for this saturation limit. 
An important consequence of the existence of e cr it is the detectability of a significant residual 
eccentricity eusA by the proposed gravitational wave detector LISA. It is found that at the 
moment of entering the LISA frequency domain eusA ~ 10~ 3 — 10~ 2 ; a signature of its 
earlier coupling with the massive circumbinary disc. We also observe large periodic inflows 
across the gap, occurring on the binary and disc dynamical time scales rather than on the 
viscous time. These periodic changes in the accretion rate (with amplitudes up to ~ 100%, 
depending on the binary eccentricity) can be considered a fingerprint of eccentric sub-parsec 
binaries migrating inside a circumbinary disc. 
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1 INTRODUCTION 

Supermassive black hole (BH) binaries are currently postulat ed to 
form in the aftermath of galaxy mergers (IBegelman et al.ll980f) . de- 
spite the difficulties, still pre sent, in identifying them observation- 
ally tsee lColpi & Dottill2009l . for a review). Thanks to advances in 
N-Body/hydrodynamical simulations, it has been shown that major 
mergers of gas-rich disc galaxies with central black h oles are con- 
duciv e to the formation of eccentric BH binaries (e.g. iMaver et al.l 
|2007|) . Orbiting inside the massive gaseous nuclear disc resulting 
upon collision, the two BHs continue to lose orbital energy and 
angular momentum under the large-scale action of gas-dynamical 
frictio n, and end up forming a circular Kep l erian binary, on parsec 
scales jEscala et al]|2005l : lDotti et alj|2007l . l2009h . As the gaseous 
and stellar mass content inside the BH orbit continues to decrease 
in response to the hardening of the binary, further inspiral is be- 
lieved to be controlled by the action of either three-body scatter- 
ing of individual stars and/or the interaction of the binary with 
a circumbinary gaseous dis c (e.g. iMerritt & Milosavlievicll200^ ; 
lArmitage & Nat araian 20q2). 
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The gravitational interaction of the massive BH binary with 
the gaseous disc is believed to be of foremost importance to as- 
sess not only its observability on sub-parsec scale, but its fate. 
Gravitational waves start to dominate the BH inspiral (leading 
to coalescence) only at tiny binary separations, of the order of 
a few milli-parsec for a binary of M ~ 1O 6 M0. If a vis- 
cous disc is present, Lindblad resonances can cause BH migra- 
tion down to the gravitational wave (GW) inspiral dom ain (e.g. 
iGoldreich & Tre maine 1980; P apaloizou & Pringlell 19771) . Follow- 
ing this proposal, a number of studies have modelled BH migra- 
tion in Keplerian, g e ometrically thin q-discs lllvanov et al.|[l999l: 



Gould & Rixl |200Cl: lArmitage & Nataraian! |2002| ; lHaiman et al.l 
20091 : lLodato et alj|2009h . 



Using high resolution hydrodynamical simulations, 



ICuadra et al.l d2009h recently investigated the evolution of the 
orbital elements of a massive BH binary, under the hypotheses (i) 
that the binary, at the radii of greatest interest (tenths of a parsec), 
is surrounded by a self-gravitating, marginally stable disc, and 
(ii) that the binary has excavated in its surroundings a cavity, i.e. 
a hollow density region of a size nearly twice the binary orbital 
separation, due to the prompt action of the binary's tidal torques. 
The simulations highlight one key aspect: that of the increase 
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of the binary eccentricity, e, during the decay of its semi-major 
axis. The excitation o f e w as already noticed and studied in 
lArmitage & Nataraianl ll2005l) . who investigated BH orbital decay 
in the presence of a Keplerian q -disc in two dimensions, as well 
as in earlier analytical work by iGoldreich & Saril d2003l) in the 
context of type-II planet migration. 

The increase of e has a number of interesting consequences. 
First, for a given semi-major axis, binaries with larger e will 
lose energy substantially faster via G Ws, coalescing on a shorter 
time scale jPeters & Mathews] 1 1963b . Second, accretion streams 
that leak through the cavity and fuel the BHs happen with a 
better defined periodicity in the case of eccentric binary (e.g., 
lArtvmowicz & Lubowlll996l) likely increasing the chance of BH 
binary identification through AGN time-variable activity. Finally, 
more eccentric binaries will retain some residual eccentricity when 
detectable by the La s er Interferometer Space Antenna (LISA ) 
teerentzenetai]|2009l : lAmaro-Seoane et al.l I20T0I : ISesanall201(jh . 
For these reasons it is important to understand if, under disc-driven 
migration, the eccentricity keeps on growing up to e ~ 1, or if there 
is a limiting eccentricity toward which the binary orbit tends. 

In this paper, we explore the binary-disc interaction with high 
resolution N-body hydr o simu lations, modelling the circumbinary 
disc as in ICuadra et all J2009h (see Section [2}- However, instead 
of starting with binaries with low eccentricities, we now construct 
a sequence of binaries with fixed semi-major axis, BH and disc- 
BH mass ratios but with different initial eccentricities eo, vary- 
ing it from 0.2 to 0.8. The binaries interact with a self-gravitating 
disc changing their orbital elements. With this approach we assess 
whether the eccentricity growth saturates, and at which value. We 
present a simple analytical interpretation of our numerical results in 
Section|4] If the saturation eccentricity is large, then the binary may 
reach coalescence with some residual eccentricity, after GW emis- 
si on has reduced it considerably . This issue was already discussed 
in lArmitage & Na taraian (2 0051) as a possible discriminant between 
gas-driven versus stellar-driven inspiral. In Section|5]we revisit this 
question in detail, in the context of the proposed LISA mission. 
The simulations also provide information on gas streams that leak 
through the cavity. We investigate how the variability properties of 
the accretion rate on to the BHs depend on the binary eccentricity. 
This analysis may lead to the identification of BH close binaries 
and estimates of their orbital elements (Section[5}- 



2 SIMULATION SET-UP 
2.1 The Model 

We model a system composed of a binary black hole surrounded 
by a gaseous disc. Since we are interested in the sub-pc separa- 
tion regime, we assume that the binary torque has already exca- 
vated an inner cavity in the gas distribution. We also assume that 
the cooling rate is long relative to the dynamical time scale, pre- 
venting disc fragmentation (e. gJRice et alj|2005h . We consider a 
binary with an initial mass ratiqj q — M2/M1 = 1/3 and a disc 
with an initial mass A/disc = 0.2M, where M — Mi + M2 is the 
total mass of the binary. The binary has initial eccentricity eo, semi- 
major axis ao, initial dynamical time td yn = f^ 1 = 2iy/Qq, where 
Slo = (GA//ao) 1//2 . Both the binary and the disc rotate in the same 
plane and direction as expected from the simulations of iDotti et al.l 



1 Unless otherwise stated, subscripts 1 and 2 refer to the primary (more 
massive) and secondary (less massive) black hole, respectively. 



j2009h . The disc is initially axisymmetric, and extends from 2ao- 
5ao. Its initial surface density profile is given by oc R~ , 

where R is the distance to the centre of mass of the system. 



2.2 Early Evolution 

ICuadra et alj d2009l) modelled the evolution of low-eccentricity bi- 
naries in the system discussed above. They found that self-gravity 
drives the initially uniformly distributed gas into a ring-like con- 
figuration located at R ~ 3ao. This ring eventually collapses and 
later spreads again in roughly the same radial range it had in the 
initial conditions (2ao-5ao). However, instead of having a uni- 
form density distrib ution, the disc displays a clear spiral pattern. 
ICuadraetai1 ( l2009l) found that this configuration remains stable for 
at least 3000^ 1 , and that during this time the binary both shrinks 
and gains eccentricity due to its interaction with the disc. In this 
study, we skip the early transient evolution and start from a snap- 
shot taken at t = 500Q.Q 1 . At this time, the disc has already settled 
into the steady-state configuration. 

2.3 The New Simulations 

Our goal is to study the secular evolution of the binary-disc sys- 
tem, focusing in the evolution of the binary eccentricity. The ideal 
method would be to follow the binary from an initial, pc-scale 
separation, until it reaches the GW-dominated regime. Unfortu- 
nately such an approach is not fe asible. The time scale for decay 
is ~ 10 4 S1 ~ 1 dCuadra et al.ll2009T) . much longer than what we can 
feasibly simulate with current computational power. Moreover, as 
the binary shrinks, its angular momentum is transferred to the disc. 
Without appropriate boundary conditions, this results in the un- 
physical expansion of the di sc, slowing further the evolution of the 
system dCuadraetai1l2009l) . To accomplish our goal we take an 
indirect approach. We run a set of simulations where the gas con- 
figuration was taken from the steady state of a previous simulation, 
as described above, but the binary had different initial eccentrici- 
ties. The energy of the binary was conserved, i.e. its semi-major 
axis a was fixed, only the angular momentum of the binary was 
changed to accomplish the various initial eccentricities eo. We then 
extrapolate the long-term evolution of the eccentricity interpreting 
the results of the different runs as snapshots of the binary life taken 
at different ages. 



2.4 Numerical method 

To simulate the binary -disc system, we use the numerical method 
described in detail by ICuadra et al.1 d2009h. We use a modified 
version of the SPH code GADGET-2 dSpringelll2005T) . We allow 
the gas to cool on a time scale which is proportional to the lo- 
cal dynamical time of the disc. To pre vent it from f r agmen ting, 
we set ft = tcooi/idyn = 10. Unlike ICuadra et all d2009h . we 
assume that the small amount of gas present in the inner cavity 
(r £j 1.75a) is isothermal, with an internal energy per unit mass 
u ~ 0.14(GAf / R). The effect of this recipe is to confine the gas 
in the inner region to a relatively thin geometry. The gravitational 
interaction between particles is calculated with a Barnes-Hut tree. 
For all runs we use 2 mi llion particles , a num ber which has been 
shown to be sufficient bv lCuadra et al.l d2009h . Since we are inter- 
ested in following the evolution of the binary orbit accurately, we 
take the BHs out of the tree and compute the gravitational forces 
acting on them directly, i.e. summing up the contributions from 
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Figure 1. Face-on view of the circumbinary disc surrounding a BH binary 
of initial eccentricity eo = 0.6 after 180 orbits. The gas density is colour- 
coded on a logarithmic scale with brighter colours corresponding to lower 
gas density; axes in units of ao- The figure shows the spiral patterns excited 
in the disc, the gap surrounding the binary, and the y in-yang sha ped gas 
inflows around the BHs. Figure made using SPLASH fpricell2007n 
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Figure 2. The eccentricity evolution of the four standard runs, starting from 
e = 0.2, 0.4, 0.6, 0.8 bottom to top. 



each gas particle. Moreover, to ensure an accurate integration, the 
dynamics of the BHs is followed with a fixed time-step, equal to 
0.0 1 f2 ( 7 1 . The BH binary is modelled as a pair of point masses, and 
their potentials are assumed to be Newtonian. Relativistic correc- 
tions, imp_ortantonly_^eji_the binary separation decays below ~ 2 
mpc dPeters & Mathewslll963h . are not included in the SPH simu- 
lations but are considered in Section|5] when estimating the eccen- 
tricity of binaries entering the LISA band. Gas particles approach- 
ing either BH are taken away from the simulation in order to avoid 
the very small time-steps they would require. They are considered 
to be accreted, and th eir mass and momentum are transfe rred to 
the corresponding BH dBate et alJl995l ; ICuadra et al.ll2006l) . In the 
present simulations the sink radius around each BH, below which 
particles are accreted, is set to 0.03ao. A face on view of the disc 
surface density is shown in Fig.[T]in which the gas has already re- 
laxed around a binary of eo = 0.6. It shows the typical spiral arms 
in the disc and the resonant streams in the inner gap region. 



3 ECCENTRICITY EVOLUTION 



As described in Section 12.31 we prepared four initial conditions 
identical but for the initial values of the binary eccentricity. In Fig. [2] 
we show the evolution of e for four runs with initial eccentricities 
eo = 0.2, 0.4, 0.6, 0.8, respectively (bottom to top). The bottom 
panel depicts the monotonic rise of e, for the run with eo = 0.2: 
The eccentricity increases almost linearly after the first 70 fg 1 . 
The run for eo = 0.4 (second panel) displays a similar behaviour, 
but the slope de/dt is much shallower (note the different scales in 
the y axes of Fig. [2}. In the third panel, corresponding to eo = 0.6, 



we observe a fast increase of the eccentricity up to e = 0.62 within 
the first few orbits; afterwards the eccentricity saturates, approach- 
ing a constant with de/dt ~ + . The top panel refers to the run 
with the largest initial eccentricity explored, eo = 0.8. This time, 
the eccentricity exhibits a negative slope with d 2 e/dt 2 steadily de- 
creasing until de/dt ~ 0~. 

The key result, illustrated in Fig. [2] is the existence of a 
limiting e cr i t that the BH binary approaches in its interaction with 
the disc. Since the runs were halted after 400 orbital cycles, we can 
only bracket the interval in which e cr it lies: e cr it € [0.62, 0.78]. 
The reason of this uncertainty is technical as we find that the 
decline of e is very hard to follow numerically due to the fast 
expulsion of the gas out of the region where torques can still 
effectively interact with the BH binary - an effect that increases 
with the binary eccentricity, as expected. Indeed, if we define i? gap 
as the inner location of the disc's half-maximal surface density, 
we find that the gas moves from an initial value of i? gap ~ 2ao 
to a time-averaged value of ~ 2.6ao, 3.0ao, 3.4ao, and 3.8ao 
during the first 53 binary orbits, for the runs with an initial binary 
eccentricity of 0.2, 0.4, 0.6, and 0.8, respectively. Such an expan- 
sion of the gas is not unexpected since no outer inflow boundary 
conditions were implemented in our simulations. While the rate of 
eccentricity change is affected by the expansion resulting from the 
initial orbital set-up, its long-term trend (whether it increases or 
decreases) is a robust conclusion from our numerical study. 

Our simulations strongly suggest the existence of a saturation 
in the disc-driven eccentricity growth, but do not pinpoint the exact 
value of e C rit- In the next section we discuss the physical reasons 
for this limit and analytically predict the value of e cr it- 
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Figure 3. Additional high eccentricity runs where the initial semi-major axis is reduced by a factor 1.8, compared to the default runs. These short runs are 
used to test the emerging picture of a limiting eccentricity depending on the amount of streams present in the cavity. The difference between the two runs is 
the fhermodynamical treatment of the gas inside the cavity. For rani this is identical to the default runs whereas in ranll we suppress gaseous inflows into the 
gap. Left panel: eccentricity versus time, for eg = 0.8. Solid line refers to rani, while dashed line to runll. Right panel: azimuthally averaged disc surface 
density as a function of R[a] in arbitrary units. Surface density is averaged over the orbits 20 — 30. Solid line refers to rani, dashed line to runll (see text for 
description). 



4 EXPLANATION OF THE SATURATION 

The growth of the eccentricity, from initial values eo below a 
critical eccentricity e cr it, and the decline of e from initial values 
eo > ficrit call for a simple physical interpretation. The increase 
of the eccentricity caused by the interaction of the binary with 
an external disc is a known fact for very unequal binaries where 
the non-axisymmetric potential pert urbations are small (as in the 
case of planetary migration, see e.g. iGoldreich & Tremaine|[l980l : 
iGoldreich & Sari 2003- |Armitage & Nataraianll2005l) . 

Goldreich & Tremain e (1980) have shown that in the high 
mass-ratio limit the binary-disc transfer of angular momentum oc- 
curs secularly through torques excited in the disc by the binary 
at discrete Lindblad and co-rotation resonances. Damping and/or 
growth of e thus depends on the relative importance of these op- 
posing torques (and so on how fluid elements are distributed in the 
disc). Principal Lindblad resonances are known to be responsible 
for opening a gap in the disc. As a consequence of disc clearance, 
co-rotation and inner Lindblad resonances are reduced in power. 
This consideration led IGoldreich & Saril (120030 to show that only 
the outer Lindblad resonances, remaining after gap opening, cause 
the increase of the eccentricity for initially low eccentric binaries. 

For the comparable mass limit studied in this paper (q = 
1/3) we have a simpler explanation. An initially small e in- 
creases because of the larger deceleration experienced by the sec- 
ondary BH near apo-a psis with respect to peri-apsis (see e.g., 
iLin & Papaloizoulll 979l ; ?). The longer time spent when nearing 
apo-apsis, and the larger over-density excited in the disc by the 
hole's gravitational pull due to its immediate proximity are both 
conducive to a net deceleration of the hole that causes the increase 
of the binary eccentricity. This increase continues as long as the 
secondary BH has a larger angular velocity at its apo-apsis cJ2, ap o 



than the fluid elements in the disc iodise- When this reverses, the 
density wake excited by the BH moves ahead imparting to the hole, 
near apo-apsis, a net tangential acceleration that tends to increase 
the angular momentum content of the binary, decreasing e. This 
argument is valid if the disc and the binary angular momenta are 
aligned. If they are antialigned (i.e. for a retrograde disc) the inter- 
action between the BHs a nd the gas increases the eccentricity up to 
e w 1 l lNixon et al.lr201 l|), resulting in a fast coalescence of the bi- 
nary. We limit our investigation to discs corotating with the binary, 
as expected if they form together duri ng a gas rich galaxy merger 
jMaver et alj|2007l ; iDotti et al.ll2009T) . In this case the torques on 
the secondary will be minimal if ajdisc = <^2, a P o- Approximating 
the binary as a purely Keplerian system and the gaseous disc to be 
in Keplerian motion around a mass Mi + Ah located at the system 
center of mass (COM), it is easy to derive: 



w 2 , a 



> 



(l + e) 2 a 3 V(l + e) 
G{A'h + Ah) 
R T 



(1) 



(2) 



where we defined Rt to be the distance of the strongest torque 
on the binary as measured from COM. Equating cjl.apo = ^disc 
yields: 

1 1 / 9 \ 

1 I , (3) 



(4) 



R T (l + e) 2 a 3 \(l + e 



which can be rearranged as 



(1 



with 8 — Rt /a. Eq. 101 implies the existence of a limiting eccen- 
tricity e cr i t that we can infer via numerical inversion of Eq. ((4). The 
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expression 

e crit = 0.66^(5-0.65) + 0.19 (5) 

provides an analytical fit to the result within a 2% accuracy in the 
range 1.8 < 5 < 4.5, relevant to our study assuming that, in a 
first approximation, 8 can be set equal to the inner edge of the disc, 
^ga P . Note that in this derivation, for a fixed S, e cr i t is indepen- 
dent of the binary mass ratio. To compare the predictions of this 
toy model with the simulations we need to define the inner edge 
of the disc. This is somewhat tricky since the disc profile is not a 
step function at a certain R/a. In our initial simulation the clean 
region within the gap has a size of i? gap ~ 2a. At larger distances 
the disc density increases reaching a maximum around R/a m 2.5. 
For 2 < S < 2.5 we get 0.55 < e crit < 0.69 which is within the 
range obtained from the numerical simulations described in Sec- 
tion [3] Note that Eq. [5] depends on the specific value of 5, i.e. on 
how close inflows of gas can get to the binary. Even though S can 
in principle be measured from our simulations, its value would also 
be affected by the lack of physical outer-boundary conditions. In- 
stead, S is usually determined equating the viscous torque in the 
accretion dis c with the positive torque exerted by the binary (see, 
e.g. Eq . 15 in lArtvmowicz & Lubow|[r994l) . lArtvmowicz & Lubowl 
dl994l) found that the size of the gap depends on e. For e « 0.6, 
q — 0.3, disc aspect ratio H/R = 0.03 and a viscous parameter 
a = 0.1, they predict i? gap ~ 2.9a, corresponding to a 5:1 com- 
mensurability resonance. Using this value for 5 we would obtain 
a larger value of e cr i t ~ 0.77. Note that the interaction between 
the binary and the disc becomes less efficient as the disc expands 
whereas the gravitational pull of the tenuous gas onto the secondary 
at peri-apsis increases. So even in a system where the influence of 
the gas inside the cavity is completely negligible, it is not clear if 
the binary could reach such a high e cr it on a relevant time scale. 
Note that a retrograde disc would not expand, since the interac- 
tion with the binary decreases its angular momentum. In this case 
the ec centricity growth remains efficient up to e « 1 jNixon et al.l 

EoTih . 

A direct com p arison between our results and 
lArtvmowicz & Lubowl dl994l) 's prediction is not straightfor- 
ward. Although our self gravitating disc is able to redistribute 
angular momentum efficiently, its total amount has to be con- 
served. Thus, discs hosting very eccentric binaries (e = 0.6, 0.8) 
keep on expanding after a short impulsive interaction with the 
binary (as discussed in Section |3}. The interaction between the 
disc and the binary is extremely inefficient when i? gap > 4a (see 
the two top panels in Fig. 2). Therefore, although a larger i? gap , 
in first approximation, implies a larger S implying a larger e cr it, it 
also results in longer timescales for the eccentricity evolution. 

The feeding of a BH binary forming in a gas rich galaxy 
merger can be a very dynamic process, and the interaction with a 
single circumbinary disc could be too idealized a picture. Larger 
scale simulations show episodic gas inflows due to the dynami- 
cal evolution of the nucleu s of the remnant (see e.g. lEscalall2006l : 
iHopkins & Ouataertll2010h . In this scenario the binary can still in- 
teract with a disc and excavate a gap, but the size of it would be 
time dependent (as in the simulations presented here) and would 
also depend on the angular momentum distribution of the inflow- 
ing streams, resulting in a range of e cr jt . 

4.1 Testing the emerging picture 

Eq. [5] shows that e cr i t depends on the location of the strongest 
torque S and thus, in first approximation, on the size of the gap 



.Rgap. In order to cross-check our results, we performed two addi- 
tional simulations of the eo = 0.8 case, in which a(l + e) was 
kept fixed, reducing the semi-major axis by a factor 1/1.8, thus in- 
creasing the relative gap size i? gap by 80%. These runs simulate a 
situation where the infalling material stays at a large distance from 
the eccentric binary and does not reach R gap ~ 2.5a, typical for 
the low eccentricity cases presented above. 

The analysis performed in the previous section only accounts 
for the pull of the disc when the secondary is at apo-apsis, neglect- 
ing torques exerted by the infalling material forming mini- accretion 
discs around the two BHs . For low eo, this approximation works 
well, because the separation of the two BHs is always much larger 
than the size of the inner mini-discs. However, in the high eo, small 
a case tested here, the secondary BH, at each peri-apsis passage, 
experiences a significant drag onto the inflowing mass accumulat- 
ing around the primary. Such drag causes the circularization of the 
orbit. Therefore, the secular evolution of the binary is determined 
by two factors: i) the distance of the gap from the secondary BH at 
apo-apsis, and ii) the amount of inflowing gas through the gap onto 
the primary BH. 

In order to separate the two effects we set two simulations with 
identical initial conditions as described above. In runl, we keep ex- 
actly the thermodynamics employed in our fiducial runs, that allow 
a stable accretion mini-disc to form around the primary hole; in 
runll the gas inside the gap evolves with the /3-cooling enabled 
just as in the rest of the disc, and can be heated by adiabatic com- 
pression. This suppresses the gaseous inflows into the cavity and 
prevents the gas from forming a significant circumprimary disc. As 
shown in the left panel of Fig. [3] rani experiences a substantial 
steady decline in e, whereas in runll, after a slight initial reduction, 
the eccentricity stays more or less constant. Such a result confirms 
our understanding of the dynamics of the system. In rani the sec- 
ondary encounters the high density region formed around the pri- 
mary at each peri-apsis passage and is slightly decelerated onto a 
more circular orbit. In runll, after a short initial relaxation phase, 
there is not enough gas in the center to cause further circularization 
(compare the two central densities in the right panel of Fig. [3}; on 
the other hand, the gap is large enough for the disc-binary interac- 
tion to be weak and, therefore, the eccentricity growth to be very 
inefficient. Note also, that for a wider binary the same effect holds, 
however, only if the secondary passes through the mini-disc of the 
primary, the size of which is independent of a. That's why in com- 
parison to the default runs in Section 3, the effect is visible more 
clearly here in the case of the narrower binary. Thus, the predicted 
limiting eccentricity e cr i t ~ 0.88 expected for 5 — 3.5 (approxi- 
mately the size of the gap in these close-separation simuations) can 
not be achieved. 

Although the torques exerted by the inside-cavity material 
when the binary eccentricity is high add complexity to the emerg- 
ing picture, this strengthen the result of a limiting eccentricity in 
the range 0.6 < e < 0.8 for the BH binary-disc configurations 
examined in this paper. 



5 OBSERVATIONAL CONSEQUENCES 

In this Section we focus on the impacts that our findings might 
have on the long-standing search for close BH binary systems in 
the Universe. First, we investigate possible periodicities residing in 
the accretion flows onto the two BHs enhancing our ability to iden- 
tify such elusive sources. Then, we study the influence of a high 
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Figure 4. Mass accretion rates onto the BHs, for the runs with eo = 
0.2,0.4,0.6,0.8 (bottom to top). Dashed (red) line refers to the primary 
BH, solid (blue) line to the lighter secondary hole. 

limiting eccentricity (attained during migration) on future gravita- 
tional wave observations with LISA. 

5.1 Periodically modulated accretion flows 

Fig. |4] shows the evolution of the accretion rate Mi and Mi onto 
each hole, for the four runs with eo = 0.2, 0.4, 0.6, 0.8. In order 
to interpret our results we consider the binary to have a total mass 
M — 3.5 x 1O 6 M0, typical for expected LISA detections, and an 
initial semi-major axis ao — 0.038 pc. Under this assumption, for a 
radiative efficiency of 0. 1 , the Eddington limit would correspond to 
accretion rates Mi, B = 0.06M Q yr" 1 and M 2 , E = 0.02M©yr _1 . 
Fig. [4] shows that this limit is fulfilled for the two high eccentric- 
ity runs only, whereas for the low-e runs the BHs accrete at super 
Eddington rates. This is possible since the numerics do not include 
any radiative feedback. The accretion rates drop significantly in the 
runs with initially higher eccentricity, owing to the expansion of the 
gap size with time (as discussed in Section 3). 

Fig. H] shows the power spectra of the accretion rates M onto 
the two BHs. The frequency / (on the x-axis) is in units of the 
binary orbital frequency fo, and the power spectral density in ar- 
bitrary units. A clear periodicity emerges at the orbital frequency 
fo, indicating a modulation of the inflow rate, induced by the or - 
bital motion dArtvmowicz & Lubowlll996l lHavasaki et al] |2008h . 
Note that smearing of the peaks at ~ fo, in Fig. [5] for eo = 0.2 
and 0.4, is due to the few-percent shrinking of the semi-major axis, 
and therefore also of the orbital period, during the evolution. As 
the binary eccentricity increases, the second and third harmonics 
of the orbital frequency increase in power and become visible. In 
the inlay of each panel the power spectrum associated to the total 
mass transfer rate onto the binary is plotted in the frequency range 
0.1/o < / < 1.0/o to illustrate the presence of other characteristic 
features at: (i) the frequency associated to the rotation of the fluid 
in the dense part of the disc: /di sc //o = (ao/f"di S c) 3//2 ; (ii) the beat 
frequency, i.e. the difference between the binary and the disc rota- 



tion frequencies: /bcat//o = l-(ao/Vdi S c) 3,/2 .Hererdi S c denotes 
the radial distance where the disc surface density has its maximum. 
Since the disc has a broad density profile, we consider the two val- 
ues r_ and r+ defined by the full width half maximum (FWHM) 
of the density and use those to estimate the expected disc and beat 
frequency intervals (enclosed by the two pairs of thick black lines 
in the inset of Fig.|5j. As expected, we observe broad features con- 
sistent with the predicted frequency ranges. Signatures of the disc 
are always visible in these plots, with a complex line structure mir- 
roring the over-densities in its spiral arms. The beat is very distinct 
in the eo = 0.8 run, marginally visible the other runs. We further 
notice that the significance of the peak in the power spectrum, at 
the binary orbital frequency fo, is weaker for low-eccentric binaries 
(eo = 0.2) than for binar ies with higher ecce ntricities. This agrees 
with previous works (cf. ICuadra et al.ll2009t) that show a mild pe- 
riodicity in the accretion rate in the case of quasi circular binaries. 
Thus, a periodic signal is expected to be a distinctive signature of 
eccentric massive BH binaries. 

The presence of periodicities in the accretion flows opens in- 
teresting prospects for monitoring sub-parsec BH eccentric binaries 
in circumbinary discs. Our fiducial system has M = 3.5 x 10 6 M© 
and an initial semi-major axis ao — 0.038 pc corresponding to an 
orbital period of 348 years, exceeding a human lifetime. Since the 
binary fingerprints in the accretion rates are related to the dynami- 
cal time, we can extrapolate our results to smaller periods as long as 
the disc and the binary are dynamically coupled (see next Section). 
For a binary with M = 3.5 x 10 6 M© and q — 1/3, binary-disc 
coupling may survive down to much shorter periods of ~ 1 month, 
making the observation of such periodicities astrophysically fea- 
sible. The interval of modulation A(M) from our runs is at the 
level of: A(M) £ [10, 50]% for e = 0.2, G [40, 100]% for 0.4, 
€ [40, 90]% for 0.6, and G [10, 50]% for 0.8. Assuming a lumi- 
nosity proportional to the time dependent accretion rate, a periodic 
monitoring of such sources will allow to construct the light curve 
for several years. An amplitude modulation of up to 100% over 
10-to-100 cycles will thus be easily identifiable. 

5.2 Residual eccentricity in GW-observations 

The existence of a limiting eccentricity that is maintained dur- 
ing the coupled evolution of the disc-binary system has impor- 
tant consequences for the detection of the binary as GW source 
in the latest stage of its evolution, i.e. during the last year of 
GW inspiral towards coalescence. Since the systems in our sim- 
ulations are far from coalescence (in our fiducial rescaling a — 
0.038pc, corresponding to ~ 10 5 Schwarzschild radii of the pri- 
mary hole), in the following we will extrapolate our findings to 
much smaller scales (order of ~ 10 3 Schwarzschild radii) making 
use of the standard opticall y thick, geometrically thin a-disc recipe 
dShakura & Sunvaevll 19731) . 

In the standard picture of BH migration, the BHs reach closer 
separations under the action of viscous torques exerted by the cir- 
cumbinary disc. This holds true as long as the migration time scale 
t m is shorter than the binary GW decay time scale tew- Since 

2 We utilize the normalized Lomb-Scargle periodogram here, wherein the 
signifi cance of each peak is directly given by the false-alarm probability 
(FAP) jScargle 1982J). Since the number of independent frequencies is the 
same for all four runs, the FAP scales identically for all runs, thus the rela- 
tive height translates into significance. For our runs, a peak needs to exceed 
a height of 12 in order to have a FAP of 0.01 
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Figure 5. Power spectrum of the accretion rate (in arbitrary units 2 ) onto the primary (blue) and secondary (red) BHs. Frequencies are in units of the initial 
binary orbital frequency fo. The inlays show zoom-ins of the power spectra, in the frequency range 0.1-l/o computed summing the accretion rate from the 
two BHs (pink). The expected intervals for the disc and the beat frequencies are marked by the thick vertical black lines, as labelled in the Figure. 



the former scales as oc a 7,/8 or q 35 ^ 16 f orgas and radiation pres- 
sure supported discs jHaiman et al .l2009h , while the latter as oc a , 
there will eventually be a critical separation a^ cc below which GW 
emission takes over and the binary decouples from the disc. Af- 
ter decoupling, binary-disc mutual torques are ineffective and the 
binary evolution is driven by GWs only. GWs tend to circularize 
the binary, but if decoupling occurs at small a, there might not be 
enough room for complete orbit circularization before entering the 
LISA frequency domain. Even a resi dual eccentricity as s mall as 
e ~ 10 - may be easily detectable JCornish & Retold) , and it 
has to be accou nted for, for a trustwor thy parameter estimation of 
the GW source jPorter & Sesanal 120101) , 

To estimate the residual eccentricity in the LISA band glisa 
we need four ingredients: 

(i) the binary eccentricity at decoupling, ed ec ; 

(ii) the binary semi-major axis at decoupling, adec', 

(iii) a model for the GW decay after decoupling; 

(iv) an estimation of /lisa at which glisa has to be computed. 

Being interested in LISA BH binaries, we consider systems charac- 
terized by 10 s M© < Mi < 10 7 M@ and 0.01 < q < 1. Item (i) is 
directly extracted from the simulations and the analytical argument 
presented in this paper. We assume that, at decoupling, the binary 
has the limiting eccentricity e cr i t given by equation (2). 

Because of its small extent, the circumbinary disc assumed in 
our simulations is unable to transfer the binary angular momen- 



tum outwards efficiently for a prolonged time scale. It is there- 
fore unsuitable for estimating a disc-driven binary decay rate to 
be compared to the GW angular momentum loss. A viable short 
cut to compute adcc (item (ii)) is to link our disc to a standard 
thin accretion disc and to estimate the gas-driven migration time 
scale in that approximation. When scaled to physical units, our bi- 
nary has ao = 0.038 pc. At such a separation, the circumbinary 
disc can be described as a stead y-state, geometrically thin, optically 
thick Shakura-Sunyaev a disc l lHaiman et alj|2009l) . Accordingly, 
the disc has a mass 

/ • \ 7/io 

M d = 1.26 x 10 3 M Q a~ o r (^) M?'*^ - R^), 

where ao.3 is viscosity parameter normalized to 0.3, rh = M/Mb 
is the accretion rate (in units of the Eddington rate), eo.i is the ra- 
diative efficiency normalized to 0.1, Mr is the total mass of the bi- 
nary in units of 10 7 Mq; the two limiting radii of the disc, Ri n and 
flout, are expressed in units of 10 3 i?sch (with Rsch = 2GM/c 2 ) 
and correspond to R ln = 2an and i? ou t = 10an, respectively. 
With this choice we infer a total disc mass Md ~ 0.25M which 
is comparable to our relaxed disc. In such a disc the time scale for 
migrat ion of the secondary BH onto the primary is given by (Eq. 
26a of lHaiman et~ai1 d2009r) ) 

t m = 1.5xl0 s y rM 7 5/8 g 8 3/8 af /16 , (7) 
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where now 03 is the binary semi-major axis in units of 10 3 i?sch 
and q B = 4g/(l + q 2 ) is the symmetric binary mass ratio. This 
time scale has to be compared with the GW decay time scale for an 
ecc entric binary which, in th e quadrupole approximation, is given 
by JPeters & Mathewslll963h 



t G w = a^- = 7.84 x 10 6 yr Mrq^atFiey 1 , (8) 
da 



where 



F(e) = (l-e 2 )- 7 / 2 ( 



1 73 2 37 4 
1+ 24 6 +96 6 



(9) 



The disc-binary decoupling occurs when taw = im, and this 
happens somewhere in the range of binary separations between 
fldcc ~ 10 2 — 10 3 i? sc h, depending on the binary mass and mass 
ratio. From that point on the dynamics of the binary is driven by 
GW emission, only. 

To address point (iii), we inte grate the Post Newtonian equa- 
tion for eccentric binaries given bv lJunker & Schaeferl(ll992l) . fol- 
lowing the eccentricity evolution down to the last stable orbit. LISA 
will be sensitive to GWs in the frequency range 10" 4 - 0.1 Hz 
and, in general, it will be able to monitor the final year of the 
binary evolution with high signal-to-noise ratio. We therefore set 
(item (iv)) /lisa = max[10 _4 Hz, /(lyr)], where /(lyr) is the 
GW frequency observed one year before the final coalescence. Note 
that the observed GW frequency is related to the rest-frame emitted 
frequency f T as / = / r /(l + z). This means that eusA, defined 
as the eccentricity of the BH binary at the time of entrance in the 
LISA band, depends on the source redshift. The 10~ 4 Hz cut-off in 
observed frequency corresponds to higher emitted frequencies as z 
increases; binaries at higher z will be caught closer to coalescence 
and will therefore show a lower residual eccentricity. 

The predicted values of eLisA, as a function of M\ for dif- 
ferent q and z, are shown in Fig. [6] Not surprisingly, the residual 
eccentricity is larger for lighter binaries (i.e., for lighter Mi) and 
smaller mass ratios q. This is simply a consequence of the scal- 
ing with M and q s of the frequency at decoupling, fdec, an d can 
be easily understood analytically as follows. By coupling the or- 
bital decay rate to the eccentricity decay rate in the quadrupole ap- 
proxim ation (sufficient for a scaling argument, IPeters & Mathews] 
Jl963h ).we get 



k 

fo 



1 + ±2i e 2 

^ T 304 c o 



870 



-3/2 



(10) 



where f T = 2/k is the frequency of the fundamental GW har- 
monic (in the rest-frame of the source) inferred from Kepler's 
law o 3 = GM/(27t/k) 2 . Eq. j 10b allows us to compute e at 
any given frequency / r , once e and f a are provided. In our case 
e = edcc ~ 0.6, and f = /dcc(adcc). If we set the value 
/lisa = 10 -4 Hz as final frequency, Eq. \W\ in the limit of small 
final e, gives 

j. 19/18 

e L iSA oc / dcc . 

The identity t m = tew requires Ode 
pling this result to Kepler's law (i.e. 

/ ~ A/f- 20/29 -33/29 
/dec oc M q s 



oc 

„3 



(11) 

m23 /29 9 22/29^ Cou _ 

oc M/7 2 ), we get 



Finally, using Eq.QT]we obtain 

7.^-0.73 -1.2 

e L isA oc M q s 



(12) 




which is basically the M and q dependence observed in Fig. [6] 
Fig. shows how this result depends on the binary eccentricity at 



Figure 6. Residual eccentricity ens A as a function of Mi, for differ- 
ent mass ratios. Each panel refers to BH binaries at different redshifts as 
labelled in the figure. In each panel, from bottom to top, curves are for 
logq = 0, -0.5, -1, -1.5, -2. 



decoupling. We see two interesting things: firstly, there is a maxi- 
mum clisa at ed C c ~ 0.4 (i.e. 6lisa is not a monotonic function of 
edcc); secondly, as long as 0.1 < edcc < 0.7, eLiSA changes only 
within a factor of ~ 2. This is a consequence of the £gw depen- 
dence on e. The higher e, the faster the GW driven evolution, and 
the larger is ad cc - Even though edcc is larger, the binary has much 
more time to circularize before entering the LISA band, showing 
a smaller residual eccentricity clisa- We note that the exact value 
of eLiSA depends on the disc properties. It is, however, interest- 
ing that a small eosA can be associated both to a fairly circular 
e c icc ~ 0.05 binary or to a binary with ed cc > 0.95. 

These results obviously depend on the assumed disc parame- 
ters. Both a lower m and a lower a would increase t m , resulting in 
a larger adcc and, in turn, in a smaller clisa- On the other hand, 
if the BHs have large spins, the radiative efficiency e may be up to 
a factor of 3 larger, acting in the opposite direction. It is however 
worth to keep in mind that tGW oc a 4 . A change of a factor of 10 
on tm will therefore result in a change of about ~ 1.8 adcc, even- 
tually influencing eLiSA only by a factor of two. We can therefore 
consider our results robust and only mildly dependent on the details 
of the disc. 



6 CONCLUSIONS 

In this paper, we explored the dynamics of sub-pc BH binaries in- 
teracting with a circumbinary gaseous disc after they have exca- 
vated a gap in the surface density distribution. We ran a sequence 
of numerical models that differ only in the initial binary eccentricity 
eo. Our aim was to study the evolution of the eccentricity in order 
to answer the following question: does the eccentricity (which is 
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Figure 7. Residual eccentricity clisa as a function of e^ cc . Red-solid 
curve refers to q = 1/3, green-dashed curve to qr = 0.1. In the figure 
the mass of the primary BH black hole is h'h = 2.6 X 10 6 M Q and the 
redshift of the binary is 2 = 1. The shaded vertical stripe brackets the 
limiting eccentricity interval found in our simulations. 



known to increase in initially circular binaries) continue to grow up 
to e — > 1 so that BH binaries in such discs reach the GW domain 
on a nearly zero angular momentum orbit, or does e saturate, and if 
so, at which value? 

The key finding is that e converges to a limiting value e cr it- 
Binaries that start with low eccentricities (eo < e cr i t ) increase 
e up to e cr i t , whereas binaries that start with high eccentricities 
(eo > e cr it) display the opposite behaviour, i.e. their eccentricity 
declines with time approaching e cr i t . Saturation rises due to the op- 
posing action of the gravitational drag experienced by the lighter, 
secondary BH in its motion near apo-apsis. For low eccentricity or- 
bits, the secondary BH excites a density wake which lags behind the 
BH at apo-apsis, causing its deceleration (and so a rise of e). The 
opposite occurs for a highly eccentric orbit: the secondary moves 
more slowly than the disc (i.e. its angular frequency is smaller than 
the angular frequency of the adjacent fluid elements) and the den- 
sity wake moves ahead of the BH path, causing a net acceleration. 
Using this simple analytical argument, the limiting eccentricity is 
independent on the binary mass ratio, but is a function of the lo- 
cation S of the inner rim of the disc from the system center of 
mass. For the range of values 2 < S < 2.5, this argument pre- 
dicts 0.55 < e cr it < 0.79, consistent with our numerical findings. 
The larger the gap size, the higher e cr it, the longer is the time scale 
on which this limit is attained. The expectation is that BH bina- 
ries, immersed in circumbinary discs, maintain a large eccentricity 
throughout the migration process. Althought in this study we have 
focused on BH binaries, the ev olution of proto-stellar binaries oc- 
curs i n a similar geometry (e.g.. lBate1l997tlArtvmowicz & Lubowl 
1 1994b and share much of the same physics. Thus, our results are 



likely relevant for the interpre tation of the observed distribution of 
binary star eccentricities (e.g. JPourbaix et al"1l2004l) . 

The existence of a relatively large limiting eccentricity in a 
BH binary, that emerges from the migration phase, has two impor- 
tant observational consequences. Firstly, the possibility of trigger- 
ing periodic inflows of gas onto the two BHs. This would enhance 
the possibility of an electromagnetic identification of a sub-parsec 
BH binary. Here we showed that periodicities occur on the dynam- 
ical time related to the Keplerian motion of the binary (depending 
on the binary parameters, from months to hundreds of years) and 
of the inner rim of the circumbinary disc, together with the beat 
frequency between the two. These features should be discernible in 
the power spectra of active nuclei, and this issue will be explored 
in detail in a forthcoming paper. Secondly, a feasible GW signature 
of a BH binary, that evolved through disc migration, is a detectable 
residual eccentricity at the time of entrance in the LISA band. In the 
case of our setup this residual e would amount to 6lisa ~ 2 x 10~ 3 
for a coalescing source at z = 1, but can be as high as clisa > 0.1 
for a lower mass, lower q binary (with M ~ 10 Mq and q < 0.1) 
at the same redshift. Thus, this study has an impact both on searches 
of periodicities in the light curves of active BHs, as well as on GW 
data stream analysis. 
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